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Abstract 

The properties of a fictitious, fermionic, many-body system based on the 
complex zeros of the Riemann zeta function are studied. The imaginary part 
of the zeros are interpreted as mean-field single-particle energies, and one 
fills them up to a Fermi energy Ep. The distribution of the total energy is 
shown to be non-Gaussian, asymmetric, and independent of Ef in the limit 
Ef — > co. The moments of the limit distribution are computed analytically. 
The autocorrelation function, the finite energy corrections, and a comparison 
with random matrix theory are also discussed. 
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According to the Riemann hypothesis, the complex zeros of the function ((s) = J2 n ~ s 
(defined for Res > 1, and by analytic continuation to the rest of the complex plane) are 
located on the critical line s M = 1/2 + iE^, with E^ real. There is a natural connexion 
between this hypothesis and quantum mechanics that originates on a spectral interpretation 
of the complex zeros s M . Indeed, a general strategy to prove the hypothesis was suggested 
by Hilbert and Polya, who proposed to find an hermitian operator whose eigenvalues are 
precisely the E^s. This operator may be viewed as the quantization of an hypothetical clas- 
sical dynamical system, which determines the "Riemann dynamics" . Although the classical 
Hamiltonian corresponding to the Riemann dynamics is not known, there are evidences for 
its existence. When the are interpreted as quantum eigenvalues, both their statistical 
properties [1-3] as well as a semiclassical interpretation of their density [4] indicate that the 
Riemann dynamics is fully chaotic and has no time-reversal symmetry. The semiclassical 
interpretation provides in fact much more detailed information (cf Eq.(2)), in particular 
concerning the classical periodic orbits. This leads to a somewhat paradoxical situation: 
although the "Riemann Hamiltonian" is not known, the detailed properties of its dynamics 
are. From the perspective of a dynamicist, the Riemann zeta thus offers the rare opportunity 
of a chaotic motion for which the relevant dynamical information is simple and explicitly 
known. Not to mention the impressive amount of existing numerical data on the complex 
zeros [2], as well as the number-theoretic background of ((s), also useful in the present 
context. 

Many aspects of the Riemann dynamics have been investigated (see [5,6] for recent review 
articles). Here we explore new facets of the problem related to a different use of the complex 
zeros as suggested by physical analogies. We consider the location E^ > of the complex 
zeros of £(s) as the single-particle levels of a fermionic many-body system. In the mean 
field approximation, the ground-state total energy is obtained by filling the single-particle 
levels from the lowest Riemann zero up to a "Fermi energy" E F . We are interested in 
the properties of such a Fermi gas. Following nuclear physics terminology, we refer to this 
"element" as the Riemannium. 

At zero temperature, the properties of the Riemannium are described by the grand 
potential 

r Ep ('Ep 

Q(E F ) = (E — E F )p(E)dE = - / J\f(E)dE . (1) 
Jo Jo 

Here p(E) = S(E — E^) is the spectral density of the complex Riemann zeros interpreted 
as quantum eigenvalues and M{E) = J p(E)dE its counting function. Q(E) corresponds to 
the sum of the imaginary part of the Riemann zeros satisfying < E^ < E F) using E F as 
reference energy [7]. It is therefore the total energy of the system. To calculate the grand 
potential we make use of the decomposition of the spectral density in smooth plus oscillatory 
parts, p = p + p. The former, as a Weyl series, has an explicit asymptotic expansion for 
large E. The oscillating term is an interference sum over the prime numbers p — 2, 3, . . . [8] 

p(E) = -Ue^jr^eMirElogp) . (2) 

The comparison of this equation with the semiclassical Gutzwiller trace formula for the 
spectral density of a dynamical system [9], p[E) = ^2 po A po cos(S po (E)/h) where po are the 
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classical periodic orbits, shows that each prime number labels an unstable periodic orbit of 
a fully chaotic system of action S p = ET p , period T p = rlogp, Lyapounov stability A p = 1, 
repetitions labeled by r, and H — 1 [4]. An unusual fact of the Riemann dynamics that plays 
an important role in what follows is the independence with respect to the energy E of the 
periods and Lyapounov exponents (notice also the minus sign in front of Eq.(2)). 

The spectral density can be integrated twice to obtain the smooth and oscillatory con- 
tributions to the grand potential, U(E F ) = -E% \og(E F /2ix) / An + 3E F /8n - 7E F /8 - 
\ogE F /A87r + cte, and 

n(E F ) - -- Re EE r i p r /nogp ■ ( 3 ) 

We are interested in the statistical properties of as a function of the Fermi energy E F . 
From Eq.(3) we have (Q) = 0. The average is done over an energy window which is small 
compared to E F but contains several oscillations of Cl (the typical scale of oscillation will be 
given below). 

As we will now see, the Riemannium has very peculiar properties. The most important 
ones demonstrated here are: i) the distribution of Cl, denoted P(Cl), is independent of 
energy as E F — > oo; ii) P(Cl) is non-Gaussian and asymmetric; iii) all the moments of 
P(Cl) may be computed with very good accuracy from Eq.(3), their value being dominated 
by the contribution of the lowest prime numbers; iv) the asymmetry of the distribution is 
due to an interference effect between repetitions of periodic orbits; v) the autocorrelation 
(Q(E F )Q(E F + e)) is a non-decaying irregular oscillatory function; vi) the finite energy 
corrections to P(Cl) are universal and well described by a circular unitary ensemble of random 
matrices (CUE). 

Some of the statistical properties of fl were explored in the past. Selberg [10] computed 
the even moments of the distribution. Up to an additive constant and a global change of 
sign, Cl(E) coincides with the function denoted Si(t) by him. More recently, Odlyzko [2] 
has numerically calculated the first four moments, but no comparison has been made with 
the analytical results of Selberg. 

On Fig. 1 is displayed the distribution P(Cl) whose properties are now discussed. We 
begin by computing the second moment. Equation (3) allows to express it as an integral 
over the period of the orbits 

(Q 2 ) = (1/2tt 2 ) / dTK(T)/T\ (4) 
Jo 

K(T) = (£ AiAj cos[E F (r t log Pl - r, log Pj )]5(T - f )>. 

Use has been made here of the definition of the form factor, K(T) = 
Ait J °° decos(eT)(p(E F )p(E F + e)), with the amplitudes defined as A { = logpj/p[ 1//2 , and 
T = (T p . + T p .)/2. Rigorous arguments valid for times T min < T < T^, where T min = log 2 
is the shortest period of the system and Th = ^rcp = log(EF/2n) is the Heisenberg time 
at Ep, and heuristic for longer ones, show [1] that the form factor of the complex zeros of 
((s) tends, as E F goes to infinity, to the corresponding function K GUF of Gaussian random 
matrices with unitary symmetry (GUE). The latter behaves as Kque(T) = T for T < Th, 
and K GUE (T) = T H for T > T H - However, the replacement K = K GUE in Eq.(4) and 
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its extrapolation to short times is of no use to understand the behavior of Cl, because the 
integral diverges. In other words, the second moment is dominated by the (non-universal) 
short periodic orbits whose contribution has to be considered explicitly. We therefore keep 
as leading-order approximation the diagonal part (i = j) in the double sum (4). This gives 
the convergent sum 

® = ^EE 4 r \ 2 " 7.9 x 10- 2 . (5) 
2ir z p J r A p r log p 

This equation exhibits two of the main properties of P(Q). The first one is its asymptotic 
independence with respect to energy. This fact is not generic for dynamical systems, since 
for example for a Fermi gas in a chaotic cavity (Q 2 ) grows linearly with the Fermi energy 
[11]. It is due to the independence of the periods and Lyapounov exponents with respect to 
energy in the Riemann dynamics. The second property is its non- universality (the sum (5) 
depends on the particular properties of the short orbits). The longer, statistically universal, 
periodic orbits provide next-to-leading order corrections to Eq.(5). For times T ^> T m j n the 
form factor Kque can be used in Eq.(4), leading to the correction (independent of the prime 
numbers) 

We have numerically checked the validity of Eq.(6) down to values of E F of the order of a 
few thousands. 

The third basic feature of the distribution is its asymmetry, as revealed by the third 
moment (fi 3 ), whose computation is now sketched. The third power of Q computed from 
Eq. (3) involves products of three cosines containing as argument the action of three different 
prime numbers. This product may be expressed as a sum of cosines involving the sum and 
differences of the actions. The term involving the sum of actions has zero average. As before, 
the remaining sum is dominated by the smallest primes. Due to the non-commensurability of 
the periods of the different primes we restrict moreover to the approximation pi = pj = pk = 
p. The oscillating factors in the sum now have the typical form cos[E(ri + rj — r k ) logp]. The 
only terms of this type having a non-zero average are those satisfying r k = ri+rj. Since there 
are three different possibilities for choosing the backward repetition r k , we finally obtain the 
convergent sum 



<^o>=-T^E E kM(n + r,)V^Tog 3 p]^ 

^ V n,rj=l 

« -5.78 x 1(T 3 . (7) 



The asymmetry of P(Cl) is therefore related to a simple but non-trivial interference phe- 
nomenon occurring between two (different or not) repetitions of a given orbit compensated 
by a backward one. The extreme case where only p = 2 and its repetitions are taken into 
account in Eq.(7) yields already 90% of its value. 

The higher moments of the distribution are computed likewise. In the limit Ep — > oo 
to leading order they are given by expressions similar to, though more complicated than, 
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Eq.(7). Table I shows a comparison of the analytical results with numerical data up to 
the sixth moment. The agreement between both calculations is astonishing. Our numerical 
results also coincide for the lowest moments with those obtained in [2]. 

An additional characterization of the Riemannium comes from the autocorrelation of the 
total energy Cn(e) = (Q(Ep) d(E F + e))/(f2o). We evaluate this quantity to leading order 
by using the expansion (3) and making a diagonal approximation. This leads to 

C (A - 1 \^ cos(erlogp) 

Due to the dominance of the shortest orbits in the sum (8) Cq is an irregularly fluctuating 
non-decaying function, as illustrated on Fig. 2. The fundamental period 27r/log2 ps 9.06 
associated to the shortest orbit can be clearly observed. 

It is interesting to make a comparison of the previous results with random matrix the- 
ory. As already emphasized, the GUE is not really appropriate because no short-time scale 
equivalent to T min is built-in in the theory (and this produces a divergence of the moments). 
We therefore concentrate on the CUE which, contrary to the Gaussian one, has an inherent 
short-time scale (although more appropriate measures, not relevant for the present discus- 
sion, have been recently proposed [12]). The following analysis reveals also some striking 
similarities existing between the Riemann zeros and eigenvalues of circular ensembles. Con- 
sider an N x N unitary matrix U describing the time-periodic dynamical evolution of a 
quantum system. We fix for simplicity the periodicity to one. The Floquet spectrum of 
U is given by the eigenvalue equation Uip a = exp(i9 a )ip a , a — 1, . . . , N. ip a are the (stro- 
boscope) eigenstates and 9 a the eigenphases. The spectral density on the unit circle is 
Pcue{@) = J2a$(@ — ®a)- The 27r-periodicity of this function leads to a decomposition in 
smooth plus oscillatory parts, pcuE = N/2n + Pcue, with 

1 oo 

Pcue(9) = -Re V TrU k exp(-ik8) . (9) 
71 k=i 

As was done for Eq.(2), a direct comparison of Eq.(9) with Gutzwiller trace formula allows 
to make a "naive" semiclassical interpretation of this equation. We look at it as a sum 
over the periodic orbits of a classical map labelled by the index k — 1, 2, . . ., having period 
T k = k, action S k = k 9, (local) energy 9, stability amplitude A k = TrU k , and H — 1 
(repetitions are degenerated with the fundamental periods). The interpretation of 9 as the 
energy satisfies, as it should, the classical relation T k = dSk/d9. We recover in this analogy 
one of the basic properties of the Riemann dynamics: the independence of the periods with 
respect to energy. Short orbits correspond to T k ks T min = 1, while long ones contributing 
to the density at the scale of the mean level spacing have period T k ~ T# = 2irp = N. 
If now an ensemble of unitary matrices is considered, Eq.(9) acquires a statistical meaning 
since the prefactors TrU k have now a distribution. It is well known through, e.g., the prime 
number theorem and the Hardy-Littlewood conjecture, that asymptotically the statistical 
properties of large prime numbers are those required to mimic random matrix statistics [1,6]. 
So for long times Eqs.(2) and (9) are very similar from a statistical point of view. This is 
not the case for times T ps T min . To have a quantitative description we look at the CUE 
distribution of the grand potential, obtained by integrating Pcue twice with respect to 9, 



5 



Q CUE = (l/7r)Re^(Trf/ fc /k 2 )exp(-ik6) . (10) 

k 

The second moment is 

(^)=(l/27r 2 )E(|Tr^| 2 )/^ 4 ■ (H) 

k 

The average is done over the ensemble of matrices. This expression is the CUE analog of 
Eq.(4) (Kcue = (\TrU k \ 2 )). Using that (\TrU k \ 2 ) = k if k < N and (\TrU k \ 2 ) = N if k > N 
[13], we get 

Vcv*)= < 3-T2^ + 0( - 1/Nl) - (12) 

to be compared to Eq.(6). As for the Riemannium, the leading term ((3)/2n 2 ps 6.1 x 
10~ 2 is controlled by short times (small values of k) in the sum (11). Because the short- 
time structure differ in both cases, the constants do not coincide. But the independence 
of the leading terms with respect to N (or energy) is not a coincidence and comes from 
the structural similarities between Eqs.(2) and (9). Moreover, the next-to-leading order 
corrections are the same if one identifies [14] the corresponding Heisenberg times N = 
log(i?i?/27r). Both corrections come from times of order Tr, where the statistical properties 
agree. This situation is in contrast with the results obtained, for instance, for the distribution 
of log £(1/2 + IE), which has been shown to agree to leading order with the (finite N) CUE 
distribution, followed by non-universal corrections [14]. 

The qualitative agreement between (Cl 2 ) and {Cl CUE ) * s ^ os ^ when considering odd powers 
of Qcue- The latter can be shown to vanish in the limit N — > oo. More generally, and 
contrary to the Riemann case, in that limit the distribution of Qcue tends to a Gaussian. 
This can be seen from Eq.(10). On the one hand we have demonstrated that flcuE is 
asymptotically dominated by the lowest contributions in the sum. But TrU k is known to 
be, for k finite and N — > oo, Gaussian independent distributed [13]. Therefore has also a 
limiting normal distribution. 

The techniques employed here to analyze the distribution of the grand potential are not 
specific to the Riemann dynamics and can be applied in a wider context to the theory of 
Fermi gases. Our results indicate that in general the grand potential of a Fermi gas, as 
well as other physical quantities derived from it, are controlled by the shortest non-universal 
classical orbits. This dominance leads, via the interference mechanisms illustrated here for 
the grand potential of the Riemannium, to non-Gaussian asymmetric distributions for the 
associated quantity. Some related results illustrating this general phenomenon were obtained 
by Tsang [15], who demonstrated non-Gaussian asymmetric distributions for the error term 
in the mean square formula of |C(l/2 + iE)\, and for the fluctuations of the number of 
lattice points inside a curve, namely a circle (see also [16]) and an hyperbola (known as the 
circle problem and Dirichlet's divisor problem, respectively). The last two problems have 
a dynamical interpretation in terms of the fluctuations of the spectral counting function of 
an associated integrable system. The mechanism leading to the asymmetry in [15] is very 
similar to the one found here. In contrast, the fluctuations of the counting function in the 
Riemann case are known to be asymptotically Gaussian [10] and are not dominated by the 
short orbits. This difference in the behavior of the counting function can be traced back to 
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the different short-time behavior of the form factor for integrable and chaotic systems, that 
produces a dominance of the short orbits in the former. We find however that integrals of the 
counting function (the grand potential being the first one, cf Eq.(l)) are always dominated 
by the short orbits, irrespective of the (chaotic or integrable) nature of the underlying 
classical dynamics. On the mathematical side, the present results can be extended in a 
straightforward manner to general Dirichlet's L-functions. 

In summary we have considered, guided by physical analogies, new properties of the Rie- 
mann zeros that reveal new aspects of their dynamical interpretation as quantum eigenvalues 
of a classically chaotic system. Our results are relevant in the theory of Fermi systems, as 
well as in the general context of quantum chaotic motion. Concerning the Riemann hy- 
pothesis and the search of the Hilbert-Polya Hamiltonian, the present study suggests that 
time-periodic dynamical evolutions have to be considered as serious candidates. 

We are specially grateful to D. Hejhal, J. Keating and A. Odlyzko for fruitful discussions 
and suggestions. The Laboratoire de Physique Theorique et Modeles Statistiques is a Unite 
Mixte de Recherche de l'Universite Paris XI and CNRS. 
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FIGURES 
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FIG. 1. Distribution of $7 computed numerically at Ep « 1.44 x 10 20 (results based on data 
from A. Odlyzko). 
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FIG. 2. The autocorrelation Eq.(8) (dotted line) compared to numerical data (dots). Ep as in 
Fig. 1. On the right part a constant eg ~ 1.29 x 10 20 should be added to the abscissa. 
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TABLES 



Moment 


Semiclassics 


Numerics 


2 


7.9290 x 1(T 2 


7.928 x 10~ 2 


3 


-5.7822 x 1(T 3 


-5.785 x 10~ 3 


4 


1.4814 x 1(T 2 


1.481 x 10~ 2 


5 


-2.7787 x 1(T 3 


-2.776 x 10~ 3 


6 


4.0007 x 10~ 3 


4.001 x 10~ 3 



TABLE I. Moments of the distribution P(£l). Numerical values are computed for the distribu- 



tion in Fig. 1. 
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